1. Field of the Invention
The present invention relates to a method of constraining by dynamic production data a fine geologic model representative of the distribution, in a heterogeneous reservoir, of a physical quantity characteristic of the subsoil structure, such as permeability or porosity.
2. Description of the Prior Art
The prior art to which reference is made hereafter is described in the following publications:
Wen, X.-H., et al.: xe2x80x9cUpscaling hydraulic conductivities in heterogeneous media: An overview. Journal of Hydrology (183)xe2x80x9d, ix-xxxii, 1996;
Renard, P.: xe2x80x9cModxc3xa9lisation des xc3xa9coulements en milieux poreux hxc3xa9txc3xa9rogxc3xa9nes: calcul des permxc3xa9abilitxc3xa9s xc3xa9quivalentesxe2x80x9d. Thxc3xa9se, Ecole des Mines de Paris, Paris, 1999;
G. de Marsily: xe2x80x9cDe I""identification des systemes hydrologiquesxe2x80x9d. Thxc3xa9se, Universitxc3xa9 Paris 6, Paris, 1976;
Hu L.-Y. et al.: xe2x80x9cConstraining a Reservoir Facies Model to Dynamic Data Using a Gradual Deformation Methodxe2x80x9d, VI European Conference on the Mathematics of Oil Recovery, Peebles, 1998;
Tarantola, A.: xe2x80x9cInverse Problem Theory: Method for Data Fitting and Model Parameter Estimationxe2x80x9d. Elsevier, Amsterdam, 1987;
Anterion F. et al.: xe2x80x9cUse of Parameter Gradients for Reservoir History Matchingxe2x80x9d, SPE 18433, Symposium on Reservoir Simulation of the Society of Petroleum Engineers, Houston, 1989;
Wen X.-H. et al.: xe2x80x9cHigh Resolution Reservoir Models Integrating Multiple-Well Production Dataxe2x80x9d, SPE 38728, Annual Technical Conference and Exhibition of the Society of Petroleum Engineers, San Antonio, 1997;
Chu L. et al. : xe2x80x9cComputation of Sensitivity Coefficients With Application to the Integration of Static and Well Test Pressure Dataxe2x80x9d, Eclipse International Forum, Milan, 1994.
Numerical simulations of flow models are widely used in the petroleum industry to develop a reservoir and to predict its dynamic behavior according to various production scenarios. The geostatistical models used to represent the geologic structure of the reservoir (permeability, porosity, etc.) require an identification of a large number of grid cells that can reach about ten millions.
To be able to carry out numerical flow simulations within reasonable computing times, common practice consists in constructing a rough simulation model by grouping together grids with different properties into macrogrids and by assigning to the macrogrids an equivalent property calculated from the local properties. This operation is referred to as upscaling.
The aim of constrained reservoir characterization is to determine the parameters of the simulation model so that the latter can reproduce the production data of the reservoir to be modelled. This parameter estimation stage is also referred to as production data fitting. The flow simulation model is thus compatible with all of the available static and dynamic data.
In common practice, the parameters of the simulation model are estimated by means of a series of trials and errors using the flow simulator.
The problem of production data fitting can also be formulated as a problem of minimizing an objective function measuring the difference between the production data observed in the field and the predictions provided by the flow simulator. Minimizing is then carried out using optimization or optimum control techniques.
A method of predicting, by means of an inversion technique, the evolution of the production of an underground reservoir, notably of a reservoir containing hydrocarbons, is for example described in U.S. Pat. No. 5,764,515 filed by the assignee.
As soon as the parameters of the simulation model are adjusted, this model can be used to simulate the present and future behavior of the reservoir. An evaluation of the in-situ reserves is thus available and a development scheme optimizing the production can be determined.
Constrained reservoir characterization thus involves multiple techniques, from geostatistical modeling to optimization problems. The introduction of the main techniques used within the scope of the xe2x80x9cinversion and upscalingxe2x80x9d coupling methodology is dealt with in the section hereafter.
Geostatistics, in its probabilistic presentation, implies that a spatial variable such as the permeability, for example, can be interpreted as a particular realization of a random function, defined by its probability law at any point in space. The increasingly common use of geostatistics by oil companies leads to the construction of fine models that can reach a large number of grid cells, In fact, geostatistics allows estimation of petrophysical properties in space from local measurements. Strictly speaking, realization of the geostatistical model has to be carried out on the scale of the measurement support, and the model thus obtained can then reach several million grid cells. Numerical flow simulation on the scale of the geostatistical model is not conceivable with the power of current computers. In order to reduce the number of grids, the grids have to be grouped together, which requires computation of the equivalent properties of the new grids as a function of the properties of the small-scale grids. This operation is referred to as upscaling.
Computation of the equivalent permeability of heterogeneous porous media has been widely studied by the community of geologists, reservoir engineers and more generally of porous media physicists.
From a mathematical point of view, the process of upscaling each directional permeability can be represented by the vectorial operator F defined by:
F:Rmxe2x86x92RMkxe2x86x92Kxe2x80x83xe2x80x83(1)
k: the permeability on the scale of the geostatistical model (dimension Rm); K:the permeability on the scale of the flow simulation model (dimension RM).
Wen et al. (1997) and Renard (1999), mentioned above, gave a review of the existing techniques from the prior art. Examples of known upscaling techniques are algebraic methods which involve simple analytical rules for plausible calculation of the equivalent permeabilities without solving a flow problem. The known method referred to as xe2x80x9cpower averagexe2x80x9d technique can be selected for example. The permeability K of block xcexa9 is equal to a power average, also called average of orderw, whose exponentw ranges between xe2x88x921 and +1:                               K          ⁢                      xe2x80x83                    ⁢                      (            w            )                          =                              (                                          1                                  mes                  ⁢                                      xe2x80x83                                    ⁢                                      (                    Ω                    )                                                              ⁢                              xe2x80x83                            ⁢                                                ∫                  Ω                                ⁢                                                      k                    w                                    ⁢                                      xe2x80x83                                    ⁢                                      ⅆ                    Ω                                                                        )                                1            /            w                                              (        2        )            
The problem of the equivalent permeability calculation thus comes down to the estimation of the exponent w allowing minimization of the error induced by upscaling (defined according to a certain criterion). For media with an isotropic log-normal distribution and a low correlation length, it is well-known that:                     w        =                  1          -                      2            α                                              (        3        )            
xcex1 being the dimension in space (xcex1=1, 2 or 3).
There are also known numerical upscaling techniques wherein calculation of the equivalent permeability involves solving the pressure p and velocity v fields of a local or global flow problem:                     {                                                                                                  -                                          k                      μ                                                        ⁢                                      xe2x80x83                                    ⁢                                      ∇                    p                                                  =                                  v                  ⁢                                      xe2x80x83                                    ⁢                  dans                  ⁢                                      xe2x80x83                                    ⁢                  Ω                                                                                                                          div                  ⁢                                      xe2x80x83                                    ⁢                                      (                    v                    )                                                  =                                  0                  ⁢                                      xe2x80x83                                    ⁢                  dans                  ⁢                                      xe2x80x83                                    ⁢                  Ω                                                                                        (        4        )            
xcexc denotes the viscosity of the flowing fluid.
The problem of geologic model updating by means of dynamic data is based on the solution of an inverse problem. This naturally poses the problem of the parameterization of the permeability field in order to allow minimization of the objective function which measures the difference (in the sense of the least squares) between the dynamic data observed in the field and the simulation results.
Parameterization of geostatistical models is a fundamental point which guarantees the success of the integration of the dynamic data into the geologic models. In fact, this integration is carried out according to an iterative procedure governed by the optimization process and disturbs an initial permeability field representative of the geostatistical model considered.
Ideally, the final permeability field must not only respect all the dynamic data taken into account in the objective function, but also must preserve the geostatistical coherence of the model (average, variogram, etc.). Observance of the dynamic data is controlled by the objective function whose value is an evaluation of the data fitting quality. Concerning the coherence of the geostatistical data, it is the parameterization of the permeability field that allows control thereof.
A known technique allowing this parameterization to be carried out is the pilot point method, which is based on the principle of the conditional geostatistical simulation applied to the Gaussian type models, described for example by de Marsily (1976) as mentioned above.
Another known technique allowing this parameterization to be carried out is the gradual deformation method. As described by Hu et al. (1998), as well as in French patents 2,780,798 and 2,795,841 and in French patent application EN-01/03,194 filed by the assignee, the gradual deformation method writes a new realization of the permeability field to be estimated, assumed to be of Gaussian type which is a linear combination of realizations independent of the random function modeling the permeability field. Permeability field k is therefore given by:                               k          ⁢                      xe2x80x83                    ⁢                      (            θ            )                          =                              ∑                          i              =              1                        n                    ⁢                      xe2x80x83                    ⁢                                    θ              i                        ⁢                          xe2x80x83                        ⁢                          k              i                                                          (        5        )            
(xcex8i)1xe2x89xa6ixe2x89xa6n: the coefficients of the linear combination, and
(ki)1xe2x89xa6ixe2x89xa6n: the independent realizations of the geostatistical model considered.
In order to preserve the geostatistical properties of the model, coefficients xcex8 must meet the normality constraint as follows:                                           ∑                          i              =              1                        n                    ⁢                      xe2x80x83                    ⁢                      θ            i            2                          =        1                            (        6        )            
Coefficients xcex8 are estimated so that the resulting permeability field k(xcex8) best reproduces the dynamic data.
Unlike the pilot point method, the gradual deformation method can be applied locally or globally. Strict observance of the geostatistical properties of the model is guaranteed by respecting the normality constraint (Equation 6) without introducing an a priori model in the objective function.
Updating a geologic model with dynamic data is based on the minimization of an objective function which measures the difference between the dynamic data observed in the field and the simulation results obtained for a set value of parameters xcex8.
Several formulations are possible to define an objective function. The formulation in the sense of the least squares is the most commonly used in the petroleum field. The objective function is thus expressed as follows:                                           J            1                    ⁢                      xe2x80x83                    ⁢                      (            θ            )                          =                              1            2                    ⁢                      xe2x80x83                    ⁢                                    (                                                d                  obs                                -                                  D                  ⁢                                      xe2x80x83                                    ⁢                                      (                    θ                    )                                                              )                        T                    ⁢                      xe2x80x83                    ⁢                      C            d                          -              1                                ⁢                      xe2x80x83                    ⁢                      (                                          d                obs                            -                              D                ⁢                                  xe2x80x83                                ⁢                                  (                  θ                  )                                                      )                                              (        7        )            
with:
dobs: the dynamic data observed in the field,
D(xcex8): the simulation results for the set value of parameters xcex8,
Cd: the covariance matrix on the observations.
As described by Tarantola (1987), a formulation better suited to the solution of improperly expressed inverse problems adds a regularization term (a priori model) in the objective function:                                           J            1                    ⁢                      xe2x80x83                    ⁢                      (            θ            )                          =                                            1              2                        ⁢                          xe2x80x83                        ⁢                                          (                                                      d                    obs                                    -                                      D                    ⁢                                          xe2x80x83                                        ⁢                                          (                      θ                      )                                                                      )                            T                        ⁢                          xe2x80x83                        ⁢                          C              d                              -                1                                      ⁢                          xe2x80x83                        ⁢                          (                                                d                  obs                                -                                  D                  ⁢                                      xe2x80x83                                    ⁢                                      (                    θ                    )                                                              )                                +                                    1              2                        ⁢                          xe2x80x83                        ⁢                                          (                                  θ                  -                                      θ                    pri                                                  )                            T                        ⁢                          xe2x80x83                        ⁢                          C              θ                              -                1                                      ⁢                          xe2x80x83                        ⁢                          (                              θ                -                                  θ                  pri                                            )                                                          (        8        )            
with:
xcex8pri: a priori estimation of parameters xcex8,
Cxcex8: the covariance matrix on the parameters.
The latter formulation of the objective function has a probabilistic interpretation. In fact, in the context of a Bayesian inversion, the a priori model is given by a probability density function of any law.
For an a priori model of Gaussian law, of average xcex8pri and of covariance Cxcex8, this probability density function is written as follows:                                           f            Θ                    ⁢                      xe2x80x83                    ⁢                      (            θ            )                          ∝                  exp          ⁢                      {                                          -                                  1                  2                                            ⁢                              xe2x80x83                            ⁢                                                (                                      θ                    -                                          θ                      pri                                                        )                                T                            ⁢                              xe2x80x83                            ⁢                              C                θ                                  -                  1                                            ⁢                              xe2x80x83                            ⁢                              (                                  θ                  -                                      θ                    pri                                                  )                                      }                                              (        9        )            
In the same context, the probability of obtaining observations dobs knowing the value of parameters xcex8, or likelihood function, can then be expressed in the following form:                                           f                          D              /              Θ                                ⁢                      xe2x80x83                    ⁢                      (                          D              =                                                d                  obs                                /                θ                                      )                          ∝                  exp          ⁢                      {                                          -                                  1                  2                                            ⁢                              xe2x80x83                            ⁢                                                (                                                            d                      obs                                        -                                          D                      ⁢                                              xe2x80x83                                            ⁢                                              (                        θ                        )                                                                              )                                T                            ⁢                              xe2x80x83                            ⁢                              C                d                                  -                  1                                            ⁢                              xe2x80x83                            ⁢                              (                                                      d                    obs                                    -                                      D                    ⁢                                          xe2x80x83                                        ⁢                                          (                      θ                      )                                                                      )                                      }                                              (        10        )            
When the flow simulation operator D is linear in relation to parameters xcex8, the a posteriori probability density function still is of Gaussian law.
Minimization of objective function J1 requires calculation of the derivatives of the simulation results in relation to the parameters to be estimated, i.e.:                                           ∂            D                                ∂            θ                          ⁢                  xe2x80x83                ⁢                  (          θ          )                                    (        11        )            
This calculation of the derivatives, essential for carrying out the minimization process under the best conditions, has been the subject of considerable work. A synthesis of the minimization process is given in the aforementioned article by Chu et al. (1994).
To date, two methods are essentially used in the petroleum industry: the numerical gradients and the gradients method. Considering its qualities in terms of numerical stability and rapidity, the gradients method has been selected for the calculation of the simulation result derivatives in relation to the parameterization of the fine geostatistical model.
The use of small letters refers to the fine geostatistical model and the use of capital letters refers to the rough simulation model. By way of example:
k denotes the permeability field on the scale of the geostatistical model, whereas K denotes the permeability field on the scale of the flow simulation model (after upscaling);
d denotes the simulation results obtained from the fine geostatistical model, whereas D denotes the simulation results obtained by means of the rough simulation model (after upscaling).
The gradients method allows calculation of the derivatives of the results of a numerical flow simulation in relation to a certain number of parameters involved in the simulation model. By way of example, it is possible to calculate the derivatives of the main production results (pressure, saturation, flow rate, etc.) in relation to the petrophysical properties (permeability, porosity, etc.) assigned to zones of the reservoir.
The gradients method is based on the derivation of the defined equations of the flow model as described by Antxc3xa9rion et al. (1989) mentioned above. These defined equations have the form of a system of non-linear equations of the following type:                     {                                                                              U                  0                                =                                  U                  ini                                                                                                                          F                  ⁢                                      xe2x80x83                                    ⁢                                      (                                          θ                      ,                                              U                        n                                            ,                                              U                                                  n                          +                          1                                                                                      )                                                  =                0                                                                        (        12        )            
xcex8: the parameters to be estimated,
Uini: initialization of the unknowns to be simulated. This initialization is calculated from the initial conditions of the partial differential equation system modelling the flow,
Un: the simulation unknowns calculated at the time tn,
Un+1: the simulation unknowns calculated at the time tn+1.
System (12) is non-linear and is generally solved by means of the Newton method based on successive linearizations of non-linear system (12) as follows:                     {                                                                              U                                      (                    0                    )                                                  =                                  U                  n                                                                                                                                                                    ∂                      F                                                              ∂                                              U                                                  n                          +                          1                                                                                                      ⁢                                      xe2x80x83                                    ⁢                                      (                                          θ                      ,                                              U                        n                                            ,                                              U                                                  (                          k                          )                                                                                      )                                    ⁢                                      xe2x80x83                                    ⁢                                      (                                                                  U                                                  (                                                      k                            +                            1                                                    )                                                                    -                                              U                                                  (                          k                          )                                                                                      )                                                  =                                                      -                    F                                    ⁢                                      xe2x80x83                                    ⁢                                      (                                          θ                      ,                                              U                        n                                            ,                                              U                                                  (                          k                          )                                                                                      )                                                                                                          (        13        )            
Calculation of the derivatives of the simulation results in relation to the parameterization xcex8 is based on the direct derivation of system (12). A new linear system, whose unknowns are the derivatives ∂Un+1/∂xcex8, results from this derivation. For each parameter xcex8i, this system is expressed in the following form:                     {                                                                                                  ∂                                          U                      0                                                                            ∂                                          θ                      i                                                                      =                                                      ∂                                          U                      ini                                                                            ∂                                          θ                      i                                                                                                                                                                                                                                                            ∂                          F                                                                          ∂                                                      U                                                          n                              +                              1                                                                                                                          ⁢                                              xe2x80x83                                            ⁢                                              (                                                  θ                          ,                                                      U                            n                                                    ,                                                      U                                                          n                              +                              1                                                                                                      )                                            ⁢                                              xe2x80x83                                            ⁢                                                                        ∂                                                      U                                                          n                              +                              1                                                                                                                                ∂                                                      θ                            i                                                                                                                +                                                                                            ∂                          F                                                                          ∂                                                      U                            n                                                                                              ⁢                                              xe2x80x83                                            ⁢                                              (                                                  θ                          ,                                                      U                            n                                                    ,                                                      U                                                          n                              +                              1                                                                                                      )                                            ⁢                                              xe2x80x83                                            ⁢                                                                        ∂                                                      U                            n                                                                                                    ∂                                                      θ                            i                                                                                                                +                                                                                            ∂                          F                                                                          ∂                                                      θ                            i                                                                                              ⁢                                              xe2x80x83                                            ⁢                                              (                                                  θ                          ,                                                      U                            n                                                    ,                                                      U                                                          n                              +                              1                                                                                                      )                                                                              =                  0                                ⁢                                  xe2x80x83                                                                                        (        14        )            
The matrix of the linear system is given by the term:                     [                                            ∂              F                                      ∂                              U                                  n                  +                  1                                                              ⁢                      xe2x80x83                    ⁢                      (                          θ              ,                              U                n                            ,                              U                                  n                  +                  1                                                      )                          ]                            (        15        )            
It is the Newton matrix of system (13) at the final iteration. The second member of this linear system is given by the term:                               [                                                    -                                                      ∂                    F                                                        ∂                                          U                      n                                                                                  ⁢                              xe2x80x83                            ⁢                              (                                  θ                  ,                                      U                    n                                    ,                                      U                                          n                      +                      1                                                                      )                            ⁢                                                ∂                                      U                    n                                                                    ∂                                      θ                    i                                                                        -                                                            ∂                  F                                                  ∂                                      θ                    i                                                              ⁢                              xe2x80x83                            ⁢                              (                                  θ                  ,                                      U                    n                                    ,                                      U                                          n                      +                      1                                                                      )                                              ]                ⁢                  xe2x80x83                                    (        16        )            
The solution of this linear system (a second member per parameter) allows obtaining all the derivatives of the simulation unknowns U in relation to the desired parameterization.
By composite derivation, it is possible to express the derivatives of the main production results D in relation to the parameterization                                           ∂            D                                ∂            θ                          =                                            ∂              D                                      ∂              U                                ⁢                      xe2x80x83                    ⁢                                    ∂              U                                      ∂              θ                                                          (        17        )            
The non-linear optimization algorithms allow calculation, according to an iterative process, of a value xcex8opt of parameters xcex8 which minimizes (locally or globally) the objective function J1 to be optimized.
The simulation results from the distribution k(xcex8opt) must allow better dynamic data fitting than those obtained from the initial distribution k(xcex8(0)). xcex8(0) denotes the value of the parameters xcex8 used to initiate the optimization process.
The objective of the iteration (k+1) of such an optimization algorithm is to determine a new estimation of parameters xcex8 according to the following principle:
xcex8(k+1)=xcex8(k)+t(k)s(k)xe2x80x83xe2x80x83(18)
Calculation of a direction: direction s(k) is the solution to a certain problem linearized at xcex8(k). The formulation of this linearized problem is based on the simulation results and on their derivatives in relation to the parameterization considered. Let:                     D        ⁢                  xe2x80x83                ⁢                  (                      θ                          (              k              )                                )                ⁢                  xe2x80x83                ⁢        and        ⁢                  xe2x80x83                ⁢                                            ∂              D                        ⁢                          xe2x80x83                        ⁢                          (                              θ                                  (                  k                  )                                            )                                            ∂            θ                                              (        19        )            
Linear seeking: interval t(k) is calculated so as to meet the descent relation:
J1(xcex8hu (k)+t(k)s(k)) less than J1(xcex8(k))xe2x80x83xe2x80x83(20)
Various optimization methods are used in the petroleum industry. Examples thereof are the deepest descent method, the Fletcher-Powell method, the Levenberg-Marquardt method and the Gauss-Newton method, which are all well-known in the art.
Updating a geologic model by dynamic data is based on the combination of various methods and techniques that are discussed above. When the geostatistical model has a reasonable size, the inversion can be carried out directly thereon without using upscaling techniques. In this context, updating is carried according to the procedure illustrated in FIG. 2.
However, when the size of the geostatistical model is too large to be used directly in the flow simulator, the use of an upscaling technique becomes mandatory. The goal of upscaling is to carry out the flow simulations on a simulation model of reduced size (referred to as rough model), thus allowing obtaining of the simulation results within a reasonable time limit. In common practice, data fitting is carried out on the rough simulation model and not on the geostatistical model. The general principle of this inversion is illustrated in FIG. 3.
Unfortunately, upon convergence of the optimization process, only the simulation model is modified and it is very difficult to return to the underlying fine geostatistical model. In fact, during the inversion process, the coherence between the initial geologic model and the simulation model is not maintained. To overcome this problem, downscaling techniques have been worked out. The aim is to determine a geologic model compatible with the constrained simulation model.
These downscaling techniques are quite substantial from a numerical point of view, in particular when the geologic model is rather large in size. They do not always allow returning from the simulation scale to the geologic scale while respecting the geostatistical constraints. Furthermore, the major drawback of these techniques is that they do not guarantee that the fine geostatistical model obtained with the downscaling technique allows respecting the dynamic data (via a flow simulation on this fine model or on a simulation model after scaling).
The many publications dealing with the problem of large geologic model fitting, including notably the aforementioned publication by Wen et al. (1997), highlight the need for a new methodology for direct updating of the fine geologic model.
The method according to the invention allows updating, by the dynamic production data, a fine geologic model representative of the distribution in the reservoir of a physical quantity characteristic of the subsoil structure (the permeability or the porosity of the reservoir rocks for example).
The method provides reservoir engineers with a methodology allowing efficient updating of geologic models as dynamic data are acquired.
The method according to the invention allows direct updating, by dynamic data, of a geologic model discretized by of a fine grid pattern representative of the distribution, in an underground reservoir, of a physical quantity characteristic of the subsoil structure: the permeability (k), the porosity ("PHgr"), etc. It comprises:
parameterization of the fine geologic model by a parameterization factor (xcex8) in order to obtain the distribution of the physical quantity in this geologic model
upscaling so as to determine the distribution of the physical quantity in a simulation model defined by a rough grid pattern;
solution, via the simulation model, of fluid flow equations to obtain simulated dynamic data ; and
determination of the analytical relations connecting variations of the simulated dynamic data and corresponding variations of parameterization factor (xcex8).
According to an implementation mode, the analytical relations connecting the variations of the simulated dynamic data and the corresponding variations of the parameterization factor of the fine geologic model are determined by combining the derivatives of the simulated dynamic data in relation to the parameterization factor on the scale of the simulation model and the derivatives of the parameterization factor of the simulation model in relation to the parameterization factor of the fine geologic model.
The simulation model is preferably first calibrated in order to reduce the error induced by upscaling, for example by carrying out the following operations:
an a priori fine geologic model representative of the model studied is selected (calibration model);
first simulation results compatible with this a priori model are directly determined;
a simulation model is determined by upscaling the fine geologic model;
second simulation results compatible with the simulation model formed, depending on upscaling parameters (c) and on simulation parameters (s), are directly determined;
calibration parameters (c, s) related to upscaling and simulation are adjusted so that the simulation results obtained from the a priori model and the simulation model are compatible.
The dynamic data are for example production data such as the pressure, the gas-oil ratio (GOR) or the fraction of water in oil.
According to an implementation mode, the parameterization parameter is selected by means of a gradual deformation or a pilot point technique.
According to an implementation mode, upscaling is carried out by means of an analytical method of a power average type or by means of a numerical method by solving a local or global flow problem.
In other words, the method according to the invention essentially includes two independent stages that can be used in an iterative process: a calibration stage and a fitting stage.
The goal of the calibration stage is to reduce the error induced by the upscaling procedure carried out to perform the flow simulation. Good calibration guarantees coherence between the fine geologic model (defined by a fine grid) and the simulation model in terms of flow. This is essential to allow reproduction of the fitting already obtained with the simulation model using the underlying fine geologic model or a rougher simulation model (modelled by a grid with larger grid cells) obtained after a new scaling operation. The calibration method of the invention is based on history matching techniques. The data to be fitted are no longer the dynamic data observed in the field, but the results of a reference simulation carried out on a given geologic model representative of the geostatistical model studied. Calibration is carried out using the simulation model obtained after scaling the reference geologic model.
The main objective of fitting is to constrain, by means of the dynamic data, directly the fine geologic model and not the simulation model. Direct parameterization of the fine geologic model is therefore performed. Upscaling is carried out on the geologic model after parameterization. Fitting involves, for example, calculation of the derivatives of the simulation results in relation to the parameterization on the scale of the fine geologic model. This allows a conventional optimization process to be used in order to directly update the fine geologic model.